統計勉強メモ

お勉強

import sys

import matplotlib.pyplot as plt
from matplotlib.pyplot import *
import numpy as np
from numpy import *
import numpy.matlib
from scipy import stats as sps

1.各種定義

標本平均・標本分散

• 高校で習う平均や分散のこと
• 今持っているデータに関して，真ん中を表す値や，散らばり具合を計算する．

母平均・母分散

• 母集団の平均や分散のこと
• 母集団：今持っているデータは母集団から適当に抜き出してきた有限の標本(サンプル)である．
• 母平均や母分散は直接知ることができない．推測するしかない．

母平均・母分散の推定に使えるものたち

母平均

• 中心極限定理
• t分布

母分散

• 不偏分散
• カイ二乗分布

計算の表記

確率変数:$x$

標本平均 :$\bar{x}=m$

標本分散 : $s^2$

特に，標本数が$n$個の場合 : $\bar{x}_n=m_n$　，$s^2_n$

母平均=期待値 : $E[x]=\mu$

母分散 : $V[x]=E[x^2]-E[x]^2=\sigma^2$

2. 基本になりそうな考え方

2.1中心極限定理

中心極限定理は母平均推定に便利なツールである． 中心極限定理の特に優れている点は，母平均の推定精度(分散)を求めることができる点．

大数の法則

ここから引用
http://mathtrain.jp/centrallimit

どんな確率分布を持つ確率変数に対しても以下が成り立つ． $$\displaystyle \lim_{ n \to \infty } \bar{x}_n = \mu $$

意味：標本平均$\bar{x}_n$は標本数$n$を大きくすると，母平均$\mu$に近づいていく．

大数の法則では，母平均$\mu$を求めたいとき標本数$n$は多ければ多いほど良い，ということが分かる．しかし，具体的に標本数$n$をいくつにすれば十分なのか？が分からない．

中心極限定理

同じくここから引用
http://mathtrain.jp/centrallimit

標本数$n$をいくつにすれば，母平均の推定精度(=標本平均の母分散)が目標に達するか？を考えたいとき，中心極限定理が登場する．

どんな確率分布(※)を持つ確率変数$x$に対しても，標本数$n$が大きいとき，以下が成り立つ． $$\bar{x}_n \sim N(\mu,\frac{\sigma^2}{n})$$

意味：標本平均の母平均$E[\bar{x}_n]$は正規分布になる．また，その正規分布の期待値は$\mu$，母分散は$\frac{\sigma^2}{n}$である．

※実は例外があるらしいby wikipedia

さらに変形すれば以下が成り立つ．

どんな確率分布を持つ確率変数$x$に対しても以下が成り立つ． $$\bar{x}_n -\mu \sim N(0,\frac{\sigma^2}{n})$$

意味：標本平均と母平均の差はゼロになることが期待される．でも実際にはゼロからばらつく．そのばらつき度合い(標本平均と母平均の差の母分散)は$\frac{\sigma^2}{n}$で与えられる．

検証

一様分布$[-\sqrt{3},\sqrt{3}]$にしたがう母集団$X$を考える．
この母集団 Xの母平均は0,母分散は1である．

1. 母集団$X$から標本数$n=1,2,\cdots,100$のそれぞれの場合についてそれぞれ標本抽出を行う．
2. それぞれの抽出結果に対して，標本平均$\bar{x}_n$を算出する(40種類の標本平均が出来上がる)
3. こうして得られる40種類の標本平均を5000セット用意する・・・①

(A)

1. 5000セットのデータから各$n$に対してヒストグラムを作成する．
2. 中心極限定理よりヒストグラムは$\mu=0$，$\sigma^2=\frac{1}{n}$の正規分布に従うと考えられる．そのような正規分布の確率密度関数をヒストグラムに併記する．
3. 5000セットのデータから各$n$に対して正規確率プロットを作成する．正規確率プロットは分布の正規性を確認する方法で，分布が正規分布ならプロットは直線上に乗る
4. 標本数$n$とShapiro-Wilk検定値との関係をプロットする

図Aから確認したいこと

• 本当に正規分布になるのか？特に標本数$n$が小さいとき．
• 中心極限定理から得られる正規分布は本当にヒストグラムを説明できるような分布か

※
正規性の検定：http://qiita.com/uri/items/e656f90e9dda342c54bb

(B)

1. 5000セットの各標本平均$\bar{x}_n$に対して不偏分散を計算する．このとき，「5000セットの不偏分散の平均値」は$\bar{x}_n$の母分散に十分に収束すると仮定する(※)．その結果，このステップで各標本平均$\bar{x}_n$の母分散が得られる・・・①
2. 中心極限定理より，各$n$に対して標本平均の母分散推定値$\sigma^2=\frac{1}{n}$を計算する・・・②
3. 横軸$n$,縦軸に分散をとり，①で得た母分散実測値と，②で得た母分散推定値をプロットする．

図Bから確認したいこと

• どんな$n$に対しても，標本平均$\bar{x}_n$の母分散が，中心極限定理によって推定可能であること．

※不偏分散
標本数が多いとき，不偏分散は母分散に収束する(大数の法則の分散バージョン)
詳しくは母分散推定の章にて．

####  A
n,m=100,5000
t=np.arange(0,n)

x=np.random.rand(n,m)*sqrt(3)*2 -sqrt(3)
nrep=np.matlib.repmat(np.arange(1,n+1),m,1).T
s=np.cumsum(x,axis=0) # sum
y=s/nrep # average

ts=np.linspace(-3,3,100)
drnum = [0,1,4,19]

coefs=[]
for i in arange(0,n):
    coefs.append(sps.shapiro(y[i,:])[1])
    if i in drnum:    
        figure()
        plt.subplot(1,2,1)
        g=plt.hist(y[i,:],label='n='+str(i+1),bins=50)
        nd=sps.norm.pdf(ts,loc=0,scale=sqrt(1/(i+1)))
        nd=nd/max(nd)*max(g[0])
        plot(ts,nd,'r')
        xlim([-3,3])
        legend()

        plt.subplot(1,2,2)
        res=sps.probplot(y[i,:],plot=plt)
        grid()

figure()
plot(t,coefs,'.-',label="shapiro")
plot([0,n],[0.05,0.05],'r-',label="thresh")
grid()
xlabel("n")
ylabel("p")
legend()
plt.show()

### B
yvar=np.var(y,axis=1,ddof=1)
figure()
plot((t+1),1/(t+1),'b',label='population',lw=3)
plot((t+1),yvar,'r.-',label='sampling')
xlabel('n')
ylabel('variance')
legend()
plt.show()

結果

(A)

• $n=1,2$あたりは正規分布には見えない．
• $n=5$を超えたあたりから十分に正規分布とみなしても良さそう
• $n$を増やせば増やすほど正規分布に近づくというわけではない

(B)

• 標本平均の中心極限定理による母分散推定値と母分散実測値は一致感が強い(一致の定義によるけれど)から中心極限定理の妥当性が確認できる．
• 中心極限定理による母分散推定の期待値は標本数$n$に依らず母分散になる．
• 中心極限定理による母分散推定の精度は標本数$n$が大きくなるほど良好になる．

検証

(C)

正規分布$N(0,5)$にしたがう母集団$X$を考える．

1. 母集団$X$から標本数$n=1,2,\cdots,100$のそれぞれの場合についてそれぞれ標本抽出を行う．
2. それぞれの抽出結果に対して，標本平均$\bar{x}_n$を算出する(40種類の標本平均が出来上がる)
3. こうして得られる40種類の標本平均を5000セット用意する
4. 5000セットのデータから各$n$に対してヒストグラムを作成する．
5. 中心極限定理よりヒストグラムは$\mu=0$，$\sigma^2=\frac{1}{n}$の正規分布に従うと考えられる．そのような正規分布の確率密度関数をヒストグラムに併記する．
6. 5000セットのデータから各$n$に対して正規確率プロットを作成する．正規確率プロットは分布の正規性を確認する方法で，分布が正規分布ならプロットは直線上に乗る
7. 標本数$n$とShapiro-Wilk検定値との関係をプロットする

図Aから確認したいこと

• 正規分布の標本平均と標本数の関係

n=100
m=5000
t=np.arange(0,n)
x=np.random.randn(n,m)*sqrt(5)
nrep=np.matlib.repmat(np.arange(1,n+1),m,1).T
s=np.cumsum(x,axis=0) # sum
y=s/nrep # average

ts=np.linspace(-5,5,100)
drnum = [0,1,2]
coefs=[]
for i in arange(0,n):
    coefs.append(sps.shapiro(y[i,:])[1])
    if i in drnum:    
        figure()
        plt.subplot(1,2,1)
        g=plt.hist(y[i,:],label='n='+str(i+1),bins=50)
        nd=sps.norm.pdf(ts,loc=0,scale=sqrt(5/(i+1)))
        nd=nd/max(nd)*max(g[0])
        plot(ts,nd,'r')
        xlim([-5,5])
        legend()

        plt.subplot(1,2,2)
        res=sps.probplot(y[i,:],plot=plt)
        grid()

figure()
plot(t,coefs,'.-',label="shapiro")
plot([0,n],[0.05,0.05],'r-',label="thresh")
grid()
xlabel("n")
ylabel("p")
legend()
plt.show()

結果

(C)

• 確率分布が正規分布の場合，標本数$n$によらず標本平均は正規分布に従うと考えても良さそう

まとめ

自分なりに中心極限定理を書き直すと以下のようになる．

どんな確率分布を持つ確率変数$x$に対しても，標本数$n$によらず以下が成り立つ． $$E[\bar{x}]=\mu$$ $$V[\bar{x}]=\sigma^2/n$$

確率分布が一様分布の場合以下が成り立つ $$\bar{x}_n \sim N(\mu,\sigma^2/n) \ \ \ (n>5)$$ (※ $n>5$は実験値)

確率分布が正規分布の場合以下が成り立つ $$\bar{x}_n \sim N(\mu,\sigma^2/n)$$

その他の確率分布の場合も，一様分布と同様に標本数$n$を十分大きくしないと標本平均を正規分布とみなせないと考えられる．

2.2 不偏分散

不偏分散は母分散推定に便利なツール．
http://mathtrain.jp/huhenbunsan

定義

不偏分散 : $u^2$
不偏分散$u^2$は標本分散から以下のように計算可能． $$u^2=\frac{n}{n-1}s^2$$

性質

$$E[u^2]=\sigma^2$$

意味：不偏分散を計算すれば標本から母分散を推定可能である．

追記
※$V[u^2]$について
http://www.crl.nitech.ac.jp/~ida/research/memo/SampleVariance/sample_variance.pdf

http://www.crl.nitech.ac.jp/~ida/research/memo/SampleVariance/sample_variance_typical.pdf

https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=4&ved=0ahUKEwji18PLnp3RAhVHfLwKHVatCQUQFggzMAM&url=http%3A%2F%2Fns1.shudo-u.ac.jp%2F~zhang%2F11.8.ppt&usg=AFQjCNE5oqBvwRVQXiNYgqCiAmrLD012Sg&sig2=lt_1ulJyciVQIG1xkBJKvw

URLによれば二つの説が確認できた．
1 $$V[u^2] \sim \frac{S_4}{n^2}-\frac{u^2}{n}$$ なお， $$S_4 = \displaystyle \sum_{ i=1 }^{n} (x_i - \bar{x})^4$$ 2 $$V[u^2]=\frac{2\sigma^4}{n-1}$$

ただ，恐らく$V[u^2]$を算出したところで特に大きな意味はないのだと思う．
なぜなら不偏分散の分布が分からないから．
不偏分散の分布が分からないなら，結局求めた不偏分散の信頼区間を知ることはできない．
できるとすれば，不偏分散の分布が正規分布になると仮定して，母分散$\sigma^2$が$u^2\pm\sqrt{V[u^2]}$の範囲内に1σ=68%の確率で存在すると主張することくらい． 詳しくはカイ二乗分布の項にて

検証

一様分布$[-\sqrt{3},\sqrt{3}]$に従う母集団$X$を考える．
この母集団$X$の母平均は0,母分散は1である．

(A)

1. 母集団$X$から標本数$n=1,2,\cdots,50$のそれぞれの場合についてそれぞれ標本抽出を行う．
2. それぞれの抽出結果に対して，標本分散$s_n^2$と不偏分散$u_n^2$を計算する．
3. 上記の計算を5000セット用意する．
4. それぞれの標本分散$s_n^2$と不偏分散$u_n^2$に対して，5000セットの標本平均を計算する．このとき，5000セットの平均は大数の法則より十分にそれぞれの期待値に収束すると仮定する．
5. 標本数$n$に対する不偏分散・標本分散の期待値をプロットする．

図から確認したいこと

• 不偏分散の期待値$E[u_n^2]$は母分散$\sigma$に等しくなること
• 標本分散の期待値$E[s_n^2]$が母分散$\sigma$に等しくないこと

(B)
不偏分散の母分散の理論式は$V[u^2]=\frac{2\sigma^4}{n-1}$で与えられる．

1. それぞれの不偏分散$u_n^2$に対して，5000セットの不偏分散を計算する．このとき，5000セットの不偏分散は不偏分散の母分散に収束すると仮定する(図Aから分かる)
2. 標本数$n$に対する不偏分散の母分散をプロットする．
3. 不偏分散の母分散の理論値をプロットする．

図Bから確認したいこと

• $V[u^2]$の傾向
• $V[u^2]$の理論式の確認

n=50
m=5000
t=np.arange(0,n)

x=np.random.rand(n,m)*sqrt(3)*2 -sqrt(3)
nrep=np.matlib.repmat(np.arange(1,n+1),m,1).T

meanrep=np.matlib.repmat(np.mean(x,axis=0),n,1)
S=np.cumsum((x-meanrep)**2,axis=0)
s2=np.mean(S/nrep,axis=1)
u2mat=S/(nrep-1)
u2=np.mean(u2mat,axis=1)

figure()
plot([0,n],[1,1],'k',label="ideal")
plot(t,s2,'r.-',label="s2")
plot(t,u2,'c.-',label="u2")
legend()
xlabel("n")
ylabel("mean of variance")
grid()

u2_var=np.var(S/(nrep-1),axis=1,ddof=1)
u2_var_t=(2)/(t-1)
figure()
plot(t,u2_var,'c.-',label='experiment')
plot(t,u2_var_t,'m-',label='theory')
grid()
xlabel("n")
ylabel("variance of variance")
legend()

plt.show()

C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in true_divide
  # This is added back by InteractiveShellApp.init_path()
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:23: RuntimeWarning: divide by zero encountered in true_divide
C:\Users\acketred\anaconda3\lib\site-packages\numpy\core\_methods.py:193: RuntimeWarning: invalid value encountered in subtract
  x = asanyarray(arr - arrmean)
C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:24: RuntimeWarning: divide by zero encountered in true_divide

検証結果

※検証ミスってる可能性あり

(A)

• $E[u_n^2]=\sigma^2$は誤り．
• $\displaystyle \lim_{ n \to \infty } E[u_n^2]=\sigma^2$は正しい．
• $E[s_n^2]\neq\sigma^2$は正しい
• 標本数$n\lt20$程度の範囲では不偏分散よりも標本分散のほうが母分散に近い

(B)

• 標本数$n$が大きくなると$V[u_n^2]$は小さくなる．
• 理論式と一致しない←理論式の証明を今度追ってみる

まとめ

検証結果より，不偏分散に関する法則を自分なりに書く．

分散に関する大数の法則(仮称)

どんな確率分布を持つ確率変数に対しても以下が成り立つ $$\displaystyle \lim_{ n \to \infty } E[u_n^2] = \sigma^2 $$ 意味：不偏分散$u_n^2$は標本数$n$を大きくすると母分散$\sigma^2$に近づいていく．

その一方で，$E[u_n^2]=\sigma^2$は必ずしも成立しないことに気をつけるべき

3. 母平均推定とt分布

母平均推定は以下のパターンに分かれる．
http://www.geisya.or.jp/~mwm48961/linear_algebra/conf1.htm

標本平均の確率分布	標本数$n<30$	標本数$n>30$
母分散が既知	$N(\mu,\sigma^2/n)$	$N(\mu,\sigma^2/n)$
母分散が未知	$t[n-1]\circ u^2$	$N(\mu,u^2/n)$

全体の流れ
母分散が未知か既知か・標本数が多いか少ないかに関わらず，母平均は標本平均$\bar{x}$を使えば尤もらしく求めることができる($E[\bar{x}]=\mu$)
しかし，問題はその標本平均$\bar{x}$の精度がどれくらいか？ということである．

この問題に答えるためには標本平均の確率分布を知る必要がある．
確率分布を用いれば，◯%の確率で標本平均±誤差範囲の中に母平均が存在する，という主張ができるようになる．

上記の表には，どうやって標本平均の確率分布を求めるか？というパターンが示してある．

3.1 母分散既知・標本数が多い場合

中心極限定理より，標本平均$\bar{x}$の確率分布は$N(\mu,\frac{\sigma^2}{n})$にしたがう．
ということは，例えば$\bar{x}\pm\frac{\sigma^2}{n}$の範囲内に1σ=68%の確率で母平均$\mu$が存在すると主張できる．

3.2 母分散既知・標本数が少ない場合

標本数$n$が少ないので，中心極限定理によって標本平均$\bar{x}$の確率分布が$N(\mu,\frac{\sigma^2}{n})$に収束している保証はない．
そこで，以下のどちらかの仮定をおくことで「母分散既知・標本数が多い」場合と同様の主張が可能である．

1. 確率変数$x$の確率分布は，少ない標本数$n$に対しても，標本平均$\bar{x}$が十分に正規分布に収束するような分布である．
2. 確率変数$x$は正規分布に従う

3.3 母分散未知・標本数が多い場合

標本数が多い場合，不偏分散の期待値について$E[u^2]=\sigma^2$が成り立つ．
また，中心極限定理より，母集団がどんな確率分布にしたがっていても標本平均$\bar{x}$は正規分布になる．
したがって，「母分散既知・標本数が多い」場合と同様の主張が可能である．

ちなみに確率変数$x$が正規分布にしたがうとき，「母分散が未知・標本数が少ない」場合と同様にt分布を利用しても良い． 標本数が多いと結局t分布も正規分布とみなせるけれど，あくまで近似せずに厳密にいきたいならt分布の利用がおすすめ．

3.2 母分散が未知・標本数が少ない場合

この場合t分布と呼ばれる分布を利用する．
t分布についてはここが分かりやすい
http://www012.upp.so-net.ne.jp/doi/biostat/CT39/distribution.pdf

t分布利用のモチベーション

母分散が未知の場合，中心極限定理によって標本平均の母分散$\frac{\sigma^2}{n}$を求めることができない．
だから，標本平均$\bar{x}$が正規分布$N(\mu,\sigma^2/n)$に従うことを利用した区間推定ができない．
そこで，「母分散が未知・標本数が多い」場合のように，母分散$\sigma^2$の代わりに，不偏分散$u^2$を使いたくなる．
このとき，以下の問題が起こる．

1. 標本数$n$が少ないときに不偏分散の期待値$E[u^2]$は母分散に一致しない(不偏分散の項で確認済み)
2. 標本数$n$が少ないときに標本平均の確率分布は正規分布にならない（中心極限定理の項で確認済み）

だから，以下の式をあてにすることはできない． $$\bar{x} \sim N(\mu,u^2/n)$$

なぜなら，上記の問題により標本平均$\bar{x}$は正規分布でない(でも$n$が十分大きければ正規分布になるような)何らかの確率分布$f$をとるからだ．　　 $$\bar{x} \sim f(\mu,u^2/n)$$

すると，標本数$n$が小さい時，標本平均$\bar{x}_n$がとる確率分布$f$はどんな分布になるか？が知りたい．
この要求に対して，昔の偉い人は確率変数$x$が正規分布にしたがう場合について調査してくれた．
その結果，t分布と呼ばれるものが出来上がった．

t分布は以下のようなものである．

標本平均$\bar{x}_n$がとる確率分布$f$は自由度$n-1$のt分布を$\sqrt{\frac{u^2}{n}}$だけ左右に拡大して，$\mu$だけ横にオフセットした分布に等しい．

あるいは

自由度n-1のt分布は以下の変数$t_{n-1}$が従う確率分布のことである． $$t_{n-1}=\frac{\bar{x}-\mu}{\sqrt{\frac{u^2}{n}}}$$

このt分布を使うと標本数$n$が小さい時，標本平均$\bar{x}_n$がとる確率分布が分かる．
だから，標本平均の信頼区間が算出できる．

気になること
以下の変数$t_n$は自由度nのt分布に従わないのだろうか． $$t_{n}=\frac{\bar{x}-\mu}{\sqrt{\frac{\sigma^2}{n}}}$$ (これなら今まで説明した4パターン全てt分布で事足りる) →どうやら違うっぽい

3.5 母平均推定への応用例

応用 1

正規分布表より，片側$3\sigma$の範囲に全体の0.4987*2=0.9974=99.74%のデータが存在することが分かる． 標本数が十分多いとき，以下の主張ができる．

1
標本数$n=n_0$のとき，99.74%の確率で以下の式は正しい． $$\bar{x}_{n0}-3\sqrt{\frac{\sigma^2}{n_0}} \leq \mu \leq \bar{x}_{n0}+3\sqrt{\frac{\sigma^2}{n_0}}$$

2
99.74%の確率で$\bar{x}_n$が$\mu\pm\Delta \mu$の範囲内に収まると主張するために必要な標本数は以下で与えられる． $$\begin{eqnarray} 3\sqrt{\frac{\sigma^2}{n}} &=& \Delta \mu\\ \iff n&=& \frac{\sigma^2} { \left( \frac{\Delta \mu}{3} \right)^2} \end{eqnarray}$$

応用 2

正規分布にしたがう母集団から標本数$n=5$の標本抽出を行ったとき$t_{n-1}=\frac{\bar{x}-\mu}{\sqrt{\frac{u^2}{n}}}$は，$ \pm2.776$の範囲に95%の確率で存在する．(自由度4のt分布の95%信頼区間=2.776)

したがって，以下の命題が正しいと主張できる．
標本平均$\bar{x}$と母平均$\mu$の差は95%の確率で$\pm 2.776\sqrt{\frac{u^2}{n}}$以内である．

検証

ある母集団$X$が平均3,分散5の正規分布にしたがうとする．ただし，この事実は実際には分からないものとする．

1. 母集団$X$から標本数$n=5$となるように標本抽出を行い，標本平均・不偏分散を計算する．
2. この標本平均と不偏分散を5000セット用意する

(A)
t分布より，標本平均$\bar{x}$と母平均$\mu$の差は95%の確率で$\pm2.776\sqrt{\frac{u_m^2}{5}}$以内である．

1. 5000セットの標本平均$\bar{x}$に対して，それぞれが$3\pm2.776\sqrt{\frac{u_m^2}{5}}$の区間内に入るかどうか検定(答え合わせ)する．
2. 区間内に入る/入らない標本平均の個数の割合を算出し，円グラフを作成する．

(B)
t分布より，標本平均$\bar{x}$と母平均$\mu$の差は95%の確率で$\pm2.571\sqrt{\frac{\sigma^2}{5}}$以内でないかという予想がある．

1. 図Aと同様の円グラフを作成する．

図A・Bから確認したいこと

• 本当に95%の割合になるかどうか
• 予想は正しいかどうか

n,m=5,5000
t=np.arange(0,n)
mu,sig2=3,5
tn_1=sps.t.ppf(0.975,n-1)
tn=sps.t.ppf(0.975,n)

x=np.random.randn(n,m)*sqrt(sig2)+mu
y=np.mean(x,axis=0)
u2=np.var(x,axis=0,ddof=1)

rangetable=(mu-tn_1*sqrt(u2/n),mu+tn_1*sqrt(u2/n))
tf=np.logical_and( rangetable[0] < y , y < rangetable[1] )
ptable=[sum(tf),m-sum(tf)]
figure()
plt.pie(ptable,autopct="%1.2f%%",colors=["gold","lightpink"] \
        ,labels=["in","out"],startangle=90,counterclock=False)
axis('equal')

rangetable=(mu-tn*sqrt(sig2/n),mu+tn*sqrt(sig2/n))
tf=np.logical_and( rangetable[0] < y , y < rangetable[1] )
ptable=[sum(tf),m-sum(tf)]
figure()
plt.pie(ptable,autopct="%1.2f%%",colors=["gold","lightpink"] \
        ,labels=["in","out"],startangle=90,counterclock=False)
axis('equal')

plt.show()

結果

(A)

• よさそう

(B)

• 予想は正しくないっぽい

4. 母分散推定と$\chi^2$分布

不偏分散の母分散は正規分布になるか？
不偏分散の項では不偏分散の母分散の傾向を確認した．
しかし，$V[u^2]$を算出したところで特に大きな意味はない．
なぜなら不偏分散の分布が分からないから．

不偏分散の分布が分からないなら，結局求めた不偏分散の信頼区間を知ることはできない．
できるとすれば，不偏分散の分布が正規分布になると仮定して，母分散$\sigma^2$が$u^2\pm\sqrt{V[u^2]}$の範囲内に1σ=68%の確率で存在すると主張することくらい．

ここで気になるのは不偏分散の母分散は正規分布になるか？という点である．

検証

一様分布$[-\sqrt{3},\sqrt{3}]$に従う母集団$X$を考える．
この母集団$X$の母平均は0,母分散は1である．

1. 母集団$X$から標本数$n=1,2,\cdots,100$のそれぞれの場合についてそれぞれ標本抽出を行う．
2. それぞれの抽出結果に対して，不偏分散$u_n^2$を計算する．
3. 上記の計算を5000セット用意する．
4. 5000セットの不偏分散を用いてヒストグラムを作成する
5. 5000セットの不偏分散に対して正規確率グラフを作成する
6. 標本数$n$に対して，正規確率グラフの相関係数をプロットする

図から確認したいこと

• $u^2$の分布がどのくらい正規分布に近いか？
• 標本数$n$と$u^2$の分布の正規性の関係

n=100
m=5000
t=np.arange(1,n)

x=np.random.rand(n,m)*sqrt(3)*2 -sqrt(3)
#mu,sig2=3,5
#x=np.random.randn(n,m)*sqrt(sig2)+mu

nrep=np.matlib.repmat(np.arange(1,n+1),m,1).T
meanrep=np.matlib.repmat(np.mean(x,axis=0),n,1)
S=np.cumsum((x-meanrep)**2,axis=0)
u2mat=S/(nrep-1)

drnum = [0,1,2,4,9,29]
coefs=[]
for i in t:
    coefs.append(sps.shapiro(u2mat[i,:])[1])
    if i in drnum:    
        figure()
        plt.subplot(1,2,1)
        g=plt.hist(u2mat[i,:],label='n='+str(i+1),bins=50)
        #xlim([-1,8])
        legend()

        plt.subplot(1,2,2)
        res=sps.probplot(u2mat[i,:],plot=plt)
        grid()

figure()
plot(t,coefs,'.-',label="shapiro")
plot([0,n],[0.05,0.05],'r-',label="thresh")
grid()
xlabel("n")
ylabel("p")
legend()
plt.show()

C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in true_divide
  if sys.path[0] == '':

検証結果

• 標本数$n$が小さいときは不偏分散$u^2$の確率分布は正規分布にならない．
• しかも確率分布は非対称
• だから，不偏分散の母分散を算出しても信頼区間の議論はあまり正しくない．
• 標本数$n$が多いときは 不偏分散$u^2$の確率分布は正規分布になることが多い

カイ二乗分布作成のモチベーション

標本数$n$が小さいとき，不偏分散の確率分布は正規分布にならない
でも，確率分布を何とかして知ることで，信頼区間の算出が行いたい．

そんな欲求に対して，確率変数$x$が正規分布に従うときのみ，カイ二乗分布が答えてくれる． カイ二乗分布とは，「正規分布に従う確率変数」の不偏分散が従う確率分布(に少し手を加えた分布)のことである．

※
上記の例では，確率変数$x$は一様分布に従っていたので結局カイ二乗分布を使っても信頼区間の推定はできない．
しかし，カイ二乗分布は「確率変数に何らかの分布を仮定すれば，その確率変数の不偏分散が従う確率分布を導出することが可能かもしれない」という示唆を与えてくれる．
もしどうしても「一様分布に従う確率変数」の不偏分散の確率分布が知りたいなら，カイ二乗分布の導出過程をきちんと追っていくことで，もしかしたら導出できるかもしれない．

カイ二乗分布導出に至るまでのストーリー(の予想)

https://www.youtube.com/watch?v=FGq9wxMwvRk

標本数$n$によって，期待値$E[u^2]$の信頼区間は変動する($n$→大，信頼区間→大) 問題はその信頼区間を具体的に知りたい=確率分布が実際にどうなっているかが知りたい，ということである．
ここで，全ての$E[u^2]$と$n$の組み合わせに対して確率分布を作るのは大変である．
だから(なのかは分からないが)$u^2 / \sigma^2$について確率分布を作ることを考える(正規化)
この操作によって，$u^2 / \sigma^2$の確率分布$f$は標本数$n$のみの関数となるはずである．

$$\frac{u^2}{\sigma^2}= \frac{\displaystyle \sum \left( x_i - \bar{x} \right)^2}{(n-1)\sigma^2} \sim f \\ \iff \frac{\displaystyle \sum \left( x_i - \bar{x} \right)^2}{\sigma^2} \sim f[n]$$

そして，昔の偉い人は，標本数$n$に依存する確率分布$f[n]$をどうにかして導出した．
その確率分布は今では自由度$n-1$の$\chi^2$分布と呼ばれていて，多くの統計界隈の人に親しまれている．

$$\frac{\displaystyle \sum \left( x_i - \bar{x} \right)^2}{\sigma^2} \sim \chi^2[n-1]$$

また，不偏分散の信頼区間を知りたい，という要求にのみ特化すれば以下の式が使いやすい. $$\displaystyle \sum \left( \frac{x_i - \bar{x}}{\sigma} \right)^2 \sim \chi^2[n-1]\\ \iff \frac{(n-1) u^2 }{\sigma^2} \sim \chi^2[n-1] \\ $$

なお，偉い人は以下で定義する母分散の不偏一致推定量$\upsilon^2$が自由度$n$のカイ二乗分布に従うことも発見している．

$$\upsilon^2= \frac{\displaystyle \sum \left( x_i - \mu \right)^2}{n}$$$$\displaystyle \sum \left( \frac{x_i - \mu}{\sigma} \right)^2 \sim \chi^2[n]\\ \iff \frac{n \upsilon^2 }{\sigma^2} \sim \chi^2[n] \\ $$

このあたりがどちらも母分散の不偏一致推定量である$\upsilon^2$，$u^2$を使ったことによるアナロジーで，$\upsilon^2$，$u^2$の違いというのが「自由度」という言葉にまとめられていて数学的に美しい匂いがする．自由度についてはその内深く考えてみたい．

あと文献によっては，以下のような説明をしているものもある．

• 正規化確率変数の二乗和は自由度$n$のカイ二乗分布に従う
• 正規化しようとしても母平均が分からないので苦し紛れに標本平均を使って正規化してみた結果，その確率変数の二乗和が自由度$n-1$のカイ二乗分布に従う

$$\displaystyle \sum \left( \frac{x_i - \mu}{\sigma} \right)^2 \sim \chi^2[n]$$$$\displaystyle \sum \left( \frac{x_i - \bar{x}}{\sigma} \right)^2 \sim \chi^2[n-1]$$

n=100
m=5000
t=np.arange(1,n)
mu,sig2=3,5

x=np.random.randn(n,m)*sqrt(sig2)+mu
nrep=np.matlib.repmat(np.arange(1,n+1),m,1).T
meanrep=np.matlib.repmat(np.mean(x,axis=0),n,1)
S=np.cumsum((x-meanrep)**2,axis=0)
u2mat=S/(nrep-1)

ts=np.linspace(0,60,200)

drnum = [1,2,4,9,29]
for i in drnum:
    figure()
    g=plt.hist(u2mat[i,:],label='n='+str(i+1),bins=50)
    chi2pdf=sps.chi2.pdf(ts,df=i+1,scale=sig2/i)
    chi2pdf=chi2pdf/max(chi2pdf)*max(g[0])
    plot(ts,chi2pdf,'r')
    legend()
plt.show()

C:\Users\acketred\anaconda3\lib\site-packages\ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in true_divide
  # Remove the CWD from sys.path while we load stuff.

結果

なんとなく惜しい感じ
スケーリングの問題か？
どこかで自由度が1ずれているような気がする．

応用

標本数$n=5$のとき，不偏分散の95%信頼区間を求める． そのために母分散の95%信頼限界値を求める．

両側で95%なので，保証したいのはカイ二乗分布の2.5% ~ 97.5%の範囲に母分散が存在することである．
自由度4のカイ二乗分布の95%信頼確率は両側それぞれ0.484，11.07である．したがって，母分散の95%信頼限界は以下で与えられる. $$\frac{(n-1) u^2 }{\sigma^2} \sim \chi^2[n-1] $$ より，

$$\frac{4u^2}{11.07} < \sigma^2 < \frac{4u^2}{0.484}$$

思ってみれば実用上これで十分だけど，一応不偏分散の信頼区間を算出する． 母分散は上記の区間に95%存在するのだから，以下が主張できる．

$\frac{4u}{11.07}-u^2$ ~ $\frac{4u}{0.484}-u^2$の範囲に母分散$\sigma^2$は95%の確率で存在する．

n,m=5,5000
t=np.arange(1,n)
mu,sig2=3,5
chi2u,chi2l=sps.chi2.ppf(0.975,n-1),sps.chi2.ppf(0.025,n-1)

x=np.random.randn(n,m)*sqrt(sig2)+mu
u2=np.var(x,axis=0,ddof=1)

rangetable=((n-1)*u2/chi2u,(n-1)*u2/chi2l)
tf=np.logical_and(rangetable[0]<sig2 , sig2<rangetable[1])
ptable=[sum(tf),m-sum(tf)]
figure()
plt.pie(ptable,autopct="%1.2f%%",colors=["gold","lightpink"] \
        ,labels=["in","out"],startangle=90,counterclock=False)
axis('equal')

plt.show()